module ssmean
	use compute
	use hmdamod
	use cimod
	implicit none
	
	integer, parameter	:: nssim=20000, nboot=2000
	integer, parameter	:: nCIgrid=40
	integer				:: ic
	integer, allocatable	:: uboots(:,:), urads(:,:)
	character(len=10)	:: cname

	real				:: cvtstat
	type(tsttype)		:: sstst
	
	contains
	
	function ssdata(rl) result(val)
		integer, parameter	:: nmethods=4
		logical	:: val
		real, allocatable	:: rl(:,:)
		integer	:: l, npop, ui(n_sample)
		real	:: x(n_sample),rls(2,nmethods,nssim),tstats(nssim)
		real, allocatable	:: datavec(:)

		if(ic==0) then
			call hmdaloaddata
		endif
		cname="data"//convtos(ic)
		ic=ic+1
		if(ic>300) then
			val=.false.
			return
		endif
		val=.true.
		
		call sethmdadata(datavec,.false.)
		datavec=pack(datavec,.not. datavec/=datavec)
		npop=size(datavec)
		print *,"npop=",npop
		datavec=datavec-sum(datavec)/npop
		
!$omp parallel do private(x,ui)		
		do l=1,nssim
			call rnund(npop,ui)
			x=datavec(ui)
			rls(:,1:3,l)=getrl_std_mean(x,tstats(l))
			rls(:,4,l)=getrl_new(x)
		enddo
		rl=setfinalrls(tstats,rls)
	end function
	
	function sstwosample(rl) result(val)
		use rlse_int
		integer, parameter	:: nmethods=4
		logical	:: val
		real, allocatable	:: rl(:,:)
		integer	:: l, npop1,npop2, ui(n_sample)
		real, allocatable	:: dv1(:),dv2(:) 		
		integer				:: ncs(n_sample)
		integer	:: n1,n2
		
		real	:: x(n_sample),xo(n_sample),rls(2,nmethods,nssim),tstats(nssim),m1,m2

		if(ic==0) then
			call hmdaloaddata
		endif
		cname="data"//convtos(ic)
		ic=ic+1
		if(ic>200) then
			val=.false.
			return
		endif
		val=.true.

		call sethmdadata(dv1,.true.)
		call sethmdadata(dv2,.true.)
		npop1=size(dv1)
		npop2=size(dv2)
		dv1=dv1-sum(dv1)/npop1
		dv2=dv2-sum(dv2)/npop2
		n1=n_sample/2
		n2=n_sample-n1

!$omp parallel do private(ui,xo,x,m1,m2)		
		do l=1,nssim
			call rnund(npop1,ui)
			xo(1:n1)=dv1(ui(1:n1))
			call rnund(npop2,ui)
			xo(n1+1:)=dv2(ui(1:n2))
			m1=sum(xo(1:n_sample/2))/(n_sample/2)
			x(1:n_sample/2)=2*(xo(1:n_sample/2)-m1)
			m2=sum(xo(n_sample/2+1:n_sample))/(n_sample/2)
			x(n_sample/2+1:n_sample)=-2*(xo(n_sample/2+1:n_sample)-m2)
			x=x+m1-m2
			rls(:,1:3,l)=getrl_std_mean(x,tstats(l))
			rls(:,4,l)=getrl_new(x)
		enddo
		rl=setfinalrls(tstats,rls)
	end function
	
	function ssMC2S(rl) result(val)
		integer, parameter	:: nmethods=4
		logical	:: val
		real, allocatable	:: rl(:,:)
		integer	:: l
		real	:: x(n_sample),xo(n_sample),rls(2,nmethods,nssim),tstats(nssim),m1,m2

		ic=ic+1
		if(ic>10) then
			val=.false.
			return
		endif
		val=.true.
!$omp parallel do private(x,xo,m1,m2)		
		do l=1,nssim
			xo=getrandomx2s(ic)
			m1=sum(xo(1:n_sample/2))/(n_sample/2)
			x(1:n_sample/2)=2*(xo(1:n_sample/2)-m1)
			m2=sum(xo(n_sample/2+1:n_sample))/(n_sample/2)
			x(n_sample/2+1:n_sample)=-2*(xo(n_sample/2+1:n_sample)-m2)
			x=x+m1-m2

			rls(:,1:3,l)=getrl_std_mean(x,tstats(l))
			rls(:,4,l)=getrl_new(x)
		enddo
		rl=setfinalrls(tstats,rls)
	end function

	function getrandomx2s(ic) result(x)
		integer	:: ic
		real	:: x(n_sample),xo(n_sample)
		xo=getrandomx(ic)
		call rnnoa(x)
		x=x/sqrt(10.0)
		x(1:n_sample/2)=x(1:n_sample/2)+xo(1:n_sample/2)
	end function
		
	function ssMC(rl) result(val)
		integer, parameter	:: nmethods=4
		logical	:: val
		real, allocatable	:: rl(:,:)
		integer	:: l
		real	:: x(n_sample),rls(2,nmethods,nssim),tstats(nssim)

		ic=ic+1
		if(ic>10) then
			val=.false.
			return
		endif
		val=.true.
!$omp parallel do private(x)		
		do l=1,nssim
			x=getrandomx(ic)
			rls(:,1:3,l)=getrl_std_mean(x,tstats(l))
			rls(:,4,l)=getrl_new(x)
		enddo
		rl=setfinalrls(tstats,rls)
	end function

	function setfinalrls(tstats,rls) result(val)
		real	:: tstats(nssim),rls(:,:,:)
		real	:: val(3,size(rls,dim=2)+1)
		real	:: cv,c,cv2
		integer	:: i
		
		val(1:2,1:size(rls,dim=2))=sum(rls,dim=3)/nssim
		cv=quantile_s(abs(tstats(:)),1-level)
		val(1:2,size(rls,dim=2)+1)=[level,val(2,1)*cv/cvtstat]

		
		cv2=quantile_s(abs(tstats),1-.99*level)
		c=2*(cv2-cv)*(val(2,1)/cvtstat)/(.01*level)
		print *,"loss func c",c		
		val(3,:)=val(2,:)+c*val(1,:)

		!do i=1,size(rls,dim=2)+1
		!	val(3,i)=quantile_s(rls(2,i,:),0.5)
		!enddo
		!val(3,size(rls,dim=2)+1)=val(3,1)*cv/1.96
		do i=2,3 
			val(i,:)=val(i,:)/val(i,size(rls,dim=2)+1)
		enddo
	end function

	function getrandomx(ic) result(x)
   		use rnfdf_int
        use rnstt_int
		use rnchi_int
		use rnnbn_int
		integer	:: ic
		real	:: x(n_sample),xm(n_sample),y(n_sample),xsi
		real(4)	:: rknbin=0.1,pnbin=0.95
		real, parameter	:: ss(10)=[1., 2.1612, 3.11805, 1.73205, 1.49071, 1.87482, 3.16487,1.0,1.0,1.0]
        select case(ic)
			case(-1)
                call rnun(xm)
                xm=merge(0.0,xm/0.5,xm>0.5)
				x=xm-0.25
				cname="flat"
            case(1)
                call rnnoa(x)
				cname="normal"
            case(2)
                call rnnoa(x)
                x=exp(x)-exp(0.5)
				cname="log-norm"
            case(3)
       			call rnfdf(4.0,5.0,x)
                x=x-1.66666666666666666
				cname="F(4,5)"
            case(4)
                call rnstt(3.0,x)
				cname="t(3)"
			case(5)
                xsi=0.4
                call rnun(xm)
                xm=xm**(-xsi)-1/(1-xsi)
                x=xm
				cname="Par 0.4"
			case(6)
                call rnnoa(x)
                x=exp(x)
                call rnnoa(y)
				call rnun(xm)
				x=merge(x,y,xm<0.5)-0.5*exp(0.5)
				cname="Mix 1"
			case(7)
                call rnnoa(x)
                x=exp(x)
                call rnnoa(y)
!				y=2*y-1
				call rnun(xm)
				x=merge(5*x,y,xm<0.05)-5*0.05*exp(0.5)
				cname="Mix 2"
            case(11)
                call rnstt(4.0,x)
                call rnstt(2.0,xm)
				x=merge(x,-abs(xm),x>0)+(1/sqrt(2.0)-0.5)
				cname="mixture of student-t(2) and student-t(4)"
            case(12)
       			call rnfdf(4.0,6.0,x)
                x=x-4.0/6
				cname="F(4,6)"
            case(8)
       			call rnun(x)
                x=x-.5
				cname="uniform"
			case(9)
				call rnnoa(x)
                call rnnoa(xm)
				x=x*xm
				cname="prod norm"
			case(10)
                xsi=0.5
                call rnun(xm)
                xm=xm**(-xsi)-1/(1-xsi)
                x=xm
				cname="Par 0.5"
		end select   
		x=x/ss(ic)	
	end function
	
	subroutine init_2s
		use tin_int
		integer	:: i
		real	:: u(n_sample)
		if(allocated(uboots)) deallocate(uboots,urads)
		allocate(uboots(n_sample,nboot),urads(n_sample,nboot))
		do i=1,nboot
			call rnun(u)
			uboots(1:n_sample/2,i)=int(u(1:n_sample/2)*n_sample/2)+1
			uboots(n_sample/2+1:n_sample,i)=n_sample/2+int(u(n_sample/2+1:n_sample)*n_sample/2)+1
		enddo
		cvtstat=tin(1-level/2,real(n_sample-2))
	end subroutine

	subroutine init_std_mean
		use tin_int
		integer	:: i
		real	:: u(n_sample)
		if(allocated(uboots)) deallocate(uboots,urads)
		allocate(uboots(n_sample,nboot),urads(n_sample,nboot))
		do i=1,nboot
			call rnun(u)
			uboots(:,i)=int(u*n_sample)+1
			urads(:,i)=merge(1,-1,u<0.5)
		enddo
		cvtstat=tin(1-level/2,real(n_sample-1))
	end subroutine
	
	function getrl_std_mean(x,tstat) result(val)
		real	:: x(n_sample),val(2,3)		! tstat, asyboot, syboot
		real	:: xbar,xm(n_sample),xb(n_sample),bstats(nboot)
		real	:: tstat,s2,bcvs(2),bcv
		integer	:: i

		xbar=sum(x)/n_sample
		xm=x-xbar
		s2=sum(xm**2)/(n_sample-1)
		val(1,1)=merge(1,0,n_sample*xbar**2/s2>cvtstat**2)
		val(2,1)=2*cvtstat*sqrt(s2/n_sample)
		
        tstat=sum(x)/sqrt(sum(x**2))
		do i=1,nboot
			xb=xm(uboots(:,i))
			bstats(i)=sum(xb)/sqrt(sum(xb**2))
		enddo
		bcvs=quantile_v(bstats,[level/2,1-level/2])
		val(1,3)=merge(1,0,tstat<bcvs(1) .or. tstat>bcvs(2))
		val(2,3)=(transcv(bcvs(2))-transcv(bcvs(1)))*sqrt(s2/n_sample)

		bcv=quantile_s(abs(bstats),1-level)
		val(1,2)=merge(1,0,abs(tstat)>bcv)
		val(2,2)=2*transcv(bcv)*sqrt(s2/n_sample)

		tstat=transcv(tstat)
	end function		
		
	function transcv(cv) result(val)
		real	:: cv,val
		val=sign(sqrt(cv**2/(1-cv**2/n_sample)),cv)
	end function
	
	subroutine	setYfromx(x,Y,s)
		real		:: x(:),Y(0:2*k)
		real		:: s
		x=sort(x)
		Y(1:k)=x(n_sample:n_sample-k+1:-1)
		Y(k+1:2*k)=-x(1:k)
		Y(0)=sum(x(k+1:n_sample-k))
		s=sqrt(sum(x(k+1:n_sample-k)**2)-Y(0)**2/(n_sample-2*k))
		Y=Y/s
	end subroutine
	
	function getrl_new(x) result(val)
		real	:: x(n_sample), val(2)
		real	:: Yo(0:2*k),scale,Zsum
		real	:: ci(2)
		
		call setYfromx(x,Yo,scale)
		ci=[getCIlb(Yo(1:k),Yo(k+1:2*k),Yo(0),sstst),-getCIlb(Yo(k+1:2*k),Yo(1:k),-Yo(0),sstst)]
		if(tstconst .and. level/=0.05) then
			if(level>0.05) then
				ci(1)=max(ci(1),getCIlb(Yo(1:k),Yo(k+1:2*k),Yo(0),tst5))
				ci(2)=min(ci(2),-getCIlb(Yo(k+1:2*k),Yo(1:k),-Yo(0),tst5))
			else
				ci(1)=min(ci(1),getCIlb(Yo(1:k),Yo(k+1:2*k),Yo(0),tst5))
				ci(2)=max(ci(2),-getCIlb(Yo(k+1:2*k),Yo(1:k),-Yo(0),tst5))
			endif
		endif

		val(1)=merge(0,1,ci(1)<0 .and. ci(2)>0)
		val(2)=scale*(ci(2)-ci(1))
	end function

end module